Jahn- Teller instability in CgHg and CeH 6 revisited. 
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The benzene cation (CsHg") has a doublet (ei 9 ) ground state in hexagonal ring (Deh) geometry. 
Therefore a Jahn- Teller (JT) distortion will lower the energy. The present theoretical study yields 
a model Hiickel-type Hamiltonian that includes the JT coupling of the ei g electronic ground state 
with the two ei g vibrational modes: in-plane ring-bending and C-C bond-stretching. We obtain the 
JT couplings from density functional theory (DFT), which gives a JT energy lowering of 970 cm" 
in agreement with previous quantum chemistry calculations. We find a non-adiabatic solution for 
vibrational spectra and predict frequencies shifts of both the benzene cation and anion, and give a 
reinterpretation of the available experimental data. 

PACS numbers: 31.50.Gh; 31.15.Ew; 31.30.Gs 



I. INTRODUCTION 

A molecular system with a degenerate electronic 
ground state will spontaneously lower its symmetry to lift 
the degeneracy. This is the Jahn-Teller (JT) Q effect, 
and is common in molecules 0. The CkH6 (benzene) 
molecule has a long history of study 0, Q • The ben- 
zene cation represents a paradigm of the dynamical JT 
effect 0. With full D^h symmetry, the highest occupied 
molecular orbital (HOMO) ei s is doubly degenerate, as 
is the lowest unoccupied molecular orbital (LUMO) eiu- 
Hence both the cation and anion are JT unstable and 
prefer lower symmetry. In the lowest (linear) approxi- 
mation, the Born-Oppcnhcimer surface has a degenerate 
loop in configuration space surrounding the D 6 h ground 
state. Everywhere on this loop, the HOMO and LUMO 
states are split by a constant JT gap. The most symmet- 
ric points on this loop represent simple T>2h distortions. 
These distortions occur in three equivalent "acute" and 
three equivalent "obtuse" forms (see the inset to Fig. 2). 
Higher order corrections lift the degeneracy on the loop, 
but zero-point energy exceeds the barrier height, creating 
a dynamic rather than a static JT distortion. 

Ah initio approaches 00, including density func- 
tional theory (DFT) [ol Il0| have been used to predict 
the ground state geometry and the JT energy lowering 
in the benzene cation. Quantum chemistry methods have 
been used to analyze vibrational spectra of the benzene 
cation [Til IT2I Il3j . A high resolution zero-kinetic-energy 
(ZEKE) photoelectron study by Linder et al. [b|[ found 
values of the frequency shift of the in-plane ring-bending 
e2 g vibrational mode. From their analysis, they deduced 
the JT energy lowering to be 208 cm -1 . They also con- 
cluded that the cation has a global minimum for the acute 
D2/1 geometry, with a local minimum for the obtuse ge- 
ometry only 8 cm -1 higher in energy. Ab initio calcu- 
lations 3 suggest 4-5 times greater energy gain. The 
present study reconciles this discrepancy. We provide a 
microscopic model Hamiltonian for the 7r-electron sys- 



tem, which includes two independent ei g modes of JT 
coupling. We derive all parameters of this model from 
DFT calculations and find vibrational spectra for both 
the cation and anion using a non-adiabatic approxima- 
tion. Our results for the cation permit a reinterpreta- 
tion of the ZEKE spectra, allowing a larger JT energy 
lowering, which agrees with both our DFT and previous 
theory 



II. MODEL HAMILTONIAN 



The starting point is a Hiickel-type Hamiltonian for the 
non-bonding 7r electrons, with a first-neighbor hopping 
integral t\ = — < i|7i|i±l > between adjacent p z atomic 
orbitals on the six C atoms. In Slater-Koster notation 
[l5j . —t\ is called the (ppn) two-site integral. In the D 6 h 
point group, all six t\ integrals are the same. The 6- 
dimensional electronic Hilbert space has single-particle 
states with a^m ei ff , e2«, b\ g symmetries, which have 
energies —2b\,—t\,t\,2t\. 

Hopping matrix elements t decrease with distance be- 
tween atoms. For small atomic displacements, t de- 
pends linearly on the C-C bond length, t\(8Ri,i+\) — 
t { [ — giSRij+i. The electron-phonon coupling constant 
gi is the same as in the Su-Schriefer-Heeger model 0] 
for polyacetylene. Next we introduce second neighbor 
hopping t<z —< i — l\H\i + 1 >= i° + <?2<5o!i to include a 
dependence on the change Sat of the Cj_i-Cj-Cj-)_i bond 
angle a« from 120°. The constant term t® does not lift 
the degeneracy of the e levels. Therefore we put t® = 
to keep the Hamiltonian as simple as possible. We use 
harmonic restoring forces K±, K2 for the bond stretching 
and angle bending vibrations. The total Hamiltonian 
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FIG. 1: (a) DFT total energy for neutral benzene with a sym- 
metric ai s distortion around the equilibrium C-C separation 
(closed circles), (b) Eigenvalue difference for the two lowest 
7T states (closed circles). The C-H bond lengths are held fixed 
for both calculations. The solid curves are best fits using the 
model Hamiltonian Eq. Q with gi and K\ adjusted. 



TL = Ti c \ + Hvib has the form 
o 



n, 



[-*? + 9iSRi,i+i\i >< i + 1| 



+g2&oti\i — 1 >< i + 1| + h.c. 
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FIG. 2: (a) DFT total energies and (b) splitting of the dou- 
blet HOMO eigenvalues of neutral benzene with a T>2h distor- 
tion and fixed C-C bond lengths as shown in the inset. The 
solid curves are best fits using the model Hamiltonian Eq. Q 
to fix the parameters <?2 and K 2. 
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The coupling constants gj_ 2 and spring constants 
K12 are calculated from DFT using the program 
NRLMOL lj| with a Gaussian basis set and the 
Perdew-Burke-Ernzerhof exchange-correlation potential 
. This yields as the optimal geometry of neutral ben- 
zene in symmetry, the C-C and C-H bond lengths of 
i?o=l-398 A and 1.08 A respectively. The energy dif- 
ference between the HOMO e\ g and a^ u orbitals defines 
the hopping integral t\ = 2.72 eV. We fix the C-H bond 
length and introduce a symmetric a\ g breathing distor- 
tion, which alters the C-C bond lengths. The results, 
shown on Fig. , can be fitted to the predictions of the 
model Hamiltonian Q E = E Q + 3Ki6R 2 (for the total 
energy) and AE = t\ = t\ — gi8R for the E{ei g ) — E(a 2u ) 
eigenvalue difference. The results are g\ = 5.27 eV/A 
and K x = 47.4 eV/1 2 . 

Second neighbor integrals are found from DFT energies 
for distortions with T>2h symmetry. We fixed the C-C 
bond lengths and varied the bond angle da as shown in 
the inset to Fig. Specifically, angles 1 and 4 were 

decreased by 2a, while the remaining four angles were 
increased by a. The hydrogen atoms were fixed along the 
bisectors of the C-C-C angles. The results were fitted to 
the predictions of the Hamiltonian Q E = Eq + QK2&a 2 
for the total energy, and Ag2&a for the splitting of doubly 
degenerate e\ g HOMO states into &2 S and bz g singlets. 
The results are K 2 = 7 '.45 eV/rad 2 and g 2 = 0.91 eV/rad. 



TABLE I: The geometry and energies of the model Hamilto- 
nian and DFT solutions for the benzene cation. The breathing 
aig mode relaxation for the anion are shown in brackets. 



III. ADIABATIC SOLUTION 

For fixed atomic coordinates, the adiabatic potential 
energy surface (APES) can easily be calculated using Eq. 

For the T>2h distortion shown in the inset of Fig. 
@, the energy depends on three variables 8R1 = SR12, 
8R2 = 8R23 and Sa2 = —Sai/2 = 8a. The elec- 
tronic contribution is the sum of the occupied eigen- 
values of Eq. (Q. The energy is not analytic at the 
point (<5i?i, <5i?2, fo)=(0,0,0), but has a cusp, character- 
istic of the first-order JT splitting. Charged benzene 
lowers its energy by 8E\ + 8E 2 + SE 3 corresponding to 
atomic relaxation via three types of vibrational modes. 
The first is the symmetric A\ g distortion of D^h symme- 
try, 8R1 — 8R2 — 8Rq and 8a = 0. The second and the 
third are two E2 g vibrational modes: the ring-bending 
and C-C bond-length distortions. It is more convenient 
to measure distortions 8R' around the relaxed breathing 
positions R = Ro + SRq — 1.416 A. The distorted geom- 
etry of benzene ions and corresponding energy gains are 
shown in Table [I] 

The geometry and energies of the DFT solution (Table 
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HJ are reasonably well reproduced by the model Hamil- 
tonian. The energy difference (3.7 meV) found between 
the acute or obtuse distortions is sensitive to the choice 
of exchange-correlation potential and smaller than the 
errors of calculations 9] . The DFT energy lowering and 
carbon atomic distortions are affected by hydrogen dis- 
tortions as well, which are left out completely in Eq. QJ. 
We did constrained geometry optimization with hydro- 
gens forced to follow carbons so that C-H bond divides 
C-C-C angle by half. The minimum energy in the con- 
strained calculations was only 4 meV higher in energy 
than in the fully optimized geometry for both acute and 
obtuse cations. The main conclusion is that the C-C 
bond-length distortions contribute the most to the JT 
energy lowering and degenerate level splitting in both 
the benzene cation and anion. 



IV. VIBRATIONAL MODEL HAMILTONIAN 

JT-coupled vibrational modes change their frequencies 
on a charged molecule. The observable frequency shifts 
serve as a measure of the JT energy lowering. We solve 
the model Hamiltonian for the vibrational spectra for 
both cation and anion. For the case of the cation, the 
Hamiltonian is projected on the two e\ g HOMO states 
with symmetries b-$ g and b 2g in the D 2 h point group, 

\b 3g >= -L (2|1 > +|2 > -|3 > -2|4 > -|5 > +|6 >) 

|6 2s >=i(|2>+|3>-|5>-|6>) (2) 

The Hamiltonian takes form U\HU\ — HqI + H z a z + 
H x a x , where U\ = (|£>3 9 >, \b 2g >), and crp are the Pauli 
matrices in the {b 3g ,b 2g } subspace, 

6 

n = ]T (i?/2M + + K 2 5a 2 j2) 

i=l 

n z = y ( 5R 'i + SR s + SR 'i + 6R 6 - 26 R 2 - 26 R s) 

+ y(6ai+<5a 4 ) 
H x = -^(S^-SR's + SR'^-SR'q) 

+ -y=(8a 2 — 5at3 + Sa 5 — 5a§) (3) 
v 12 

In adiabatic approximation, the JT orbital splitting is 
2A = 2g\/ZK 1 + 2g\jZK 2 and the JT energy lowering is 
_Ejt = — A/2. To solve the problem beyond the adiabatic 
approximation we have to quantize vibrational motions. 
We rewrite the Hamiltonian J3J) in Cartesian coordinates 
(Sxi,Syi). The kinetic energy term is diagonal. Choosing 
as the unit of displacement R= 1.416 A(the benzene ring 
radius after a\ g mode relaxation), and as the energy unit 
K1R 2 , the potential energy of the Hamiltonian H (Eq. 
©) depends on a single parameter K2/K1R 2 = A. 



After diagonalizing Hq we get twelve normal modes 
whose amplitudes we designate ©i..i2- The potential 
energy is independent of the motion of the center of 
mass {X = ~^2xi,Y = Y^Vi) anc ^ the rm S rotation 
(a = ^2 a i)- This forces three normal modes to have fre- 
quencies zero. Since we neglect C-H stretching motions, 
the mass M is M c + M H . Choosing A = 0.0784 the nine 
remaining normal modes are: the degenerate ring bend- 
ing E 2g mode at uj\ = 581 cm -1 with some admixture 
of the C-C bond length alteration, the pure ring bend- 
ing B\ u mode at 962 cm -1 , the breathing A\ g mode at 
992 cm -1 , the mostly C-C bond-length-alteration dou- 
blets E\ u and E 2g at 1262 cm -1 and lo 2 — 1644 cm -1 
and the pure C-C bond-stretching mode B 2u at 1718 
cm -1 . It is interesting to compare the model Hamil- 
tonian normal mode spectrum with the DFT NRLMOL 
calculations. There are 21 in-plane vibrational modes 
r = 2A lg + A 2g + AE 2g + 2B lu + 2B 2u + ZE lu . Six of them 
are high frequency C-H modes A\ g + B lu + E lu + E 2g 
around 3100 cm -1 . From the remaining 15 vibrations 
we identify the ring-bending E 2g mode at 600 cm -1 , the 
breathing A\ g at 992 cm" 1 , and the bond-stretching E 2g 
at 1589 cm" 1 . 

The off diagonal part H x of the Hamiltonian J2J de- 
pends only on the two JT-active E 2g vibrational normal 
modes. Let us label the amplitudes of these modes by 
(8i a ,0ib) and (82a, ©2b) and denote their frequencies 
by lui (ring bending) and oj 2 (bond-stretching) as above. 
We also use the coordinates Oi, 2 to denote the to- 
tal displacements Q\ = 9f a + Qj b etc. Then Eq. © 
becomes: 

w = Y(ef + e? 6 + w?e?) 
+ y(eL + ©26 + ^62) 

H z = -Gi©i cos(/?i) - G 2 2 cos(/3 2 ) 

H x = -G 1 © 1 sin(/3 1 )-G 2 2 sin(/3 2 ), (4) 

where the angles Pi, (3 2 contain arbitrary additive con- 
stants which are the choices of the orientation angles of 
the orthogonal eigenvectors in the two-dimensional E g 
subspaces of the two modes. Note that the adiabatic elec- 
tronic eigenvalues which diagonalize Ti. z cr z +Ti. x a x (in the 
approximation that (®i a , ..,@2b) are classical variables) 
are ±Ajt where 

Aj T = G\Q\ + Gjej + 2GiG 2 0i0 2 cos(#i - /%) (5) 

This does not depend separately on the angles /3i, (3 2 . 
The modified electron-phonon coupling constants Gi , G 2 
depend on the JT couplings 171 , 172 and on the E g normal 
mode eigenstates. For our case (A = 0.0784) we find 

Gi = T101/6 + T 2 g 2 /2R = 1.3 eV/A 

G 2 = T m /6-T 4 g 2 /2R = 4A4 eV/A (6) 

Let us consider further the adiabatic eigenstates of 
H = Ti z (T z + H x <j x . What happens when we change 
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the angle 5/3i = <5tan _1 Oib/©i a ? The rotation ma- 
trix R = cos(a/2)(T 2 +isin(a/2)cj a . has the property that 
RUB) gives a new matrix Tt' with angles j3\ — > /3i + a 
and /?2 — > /?2 + «• The energies ±Ajt are unchanged 
and the eigenfunctions are rotated in {b 3g , b 2g ) space by 
a/2. This type of change of distortion pattern is called a 
"pseudo-rotation" . A full pseudo-rotation by a = 2n re- 
stores the molecular geometry and changes the sign of the 
electronic eigenstate. This sign change under adiabatic 
distortion is the Berry phase, and needs to be correctly 
incorporated when finding the system's total wavefunc- 
tion in adiabatic approximation. In particular it leads 
to fractional quanta of pseudorotational angular momen- 
tum. Our method will find total wavefunctions without 
making the adiabatic approximation and will therefore 
yield correct answers without explicit mention of Berry 
phase. 

V. NON ADIABATIC SOLUTION 




FIG. 3: Vibrational spectrum of the Jahn- Teller unstable ion 
versus parameter a, where coupling constants rei = 1.18 * a 
and K2 = 0.8 * a. For the cation, the coupling constant for 
the ring bending mode is /ci = 1.18 and coupling to the bond 
stretching mode is k,2 = 0.8. Therefore for a = 1 the neu- 
tral molecule vibrational frequency u\ = 581 cm' 1 becomes 
0.55aji and LU2 = 1644 cm -1 reduces to 0.84^2 for the cation. 



The quantized vibrational Hamiltonian Q depends on 
three parameters uj\/lo 2 — 0.353 and the two electron- 
phonon coupling constants K\ = (Ai/hwi) 1 / 2 = 1.18, 



k 2 = (A 2 /fiw- 



a/2 



0.8, where Ai 



G\IMuj\ 



817cm" 1 and A 2 = G\jMj\ = 1052cm" 1 are twice 
the JT energy gain due to the ring-bending and the 
bond stretching modes correspondingly. Similar values 
of Ki = 1.06 and k 2 = 0.817 were obtained by Eiding et 
al. |13| using ab initio many body techniques. 

To solve Eq. Q for the vibrational spectrum, assume 
for the moment a zero coupling to the bond stretching 
vibrational mode G 2 = 0. Then Eq. depends on a 
single parameter k±. The classical solution of Longuet- 
Higgins et al. p| uses polar coordinates, which simplifies 
equations. Instead we seek the solution of Eq. (0} in the 
form 



|*>=E E 

m— no—0 



(4) W1 & 



A, 




|0 > \b 3g > 



> \b 2g > 



(7) 



Note that we are not making the adiabatic approxima- 
tion where the electronic wavefunction ty c \ is strictly held 
in the lower of the two Jahn- Teller split levels. In adia- 
batic approximation the double valuedness of \l/ c i under 
rotation by (3 = 2tt is compensated by a restriction on 
the vibrational wavefunctions. This effect is automati- 
cally handled in our non- adiabatic treatment. The oper- 
ators c\, c\ create harmonic oscillator states in the two- 
dimensional manifold of e g vibrations, with their origin 
at the symmetric {D^h) benzene ion coordinates. Al- 
though different from the basis of Longuet-Higgins et al. 
, both basis sets are complete and converge rapidly to 
the same solution. The solution is found by exact di- 
agonalization of Eq. J2J for N = 10. The eigenvalues 



coincides with those of Linder et. al. 14] for the case of 
zero quadratic coupling. 

The vibrational frequencies on the neutral molecule 
u>i = 581 cm" 1 and w 2 = 1644 cm" 1 become (for the 
benzene cation) 0.54wi and 0.7^2 in the uncoupled-mode 
approximation. To take into account the interaction be- 
tween the two Ei g vibrational modes (when both G\ ^ 
and G 2 ^ 0) one has to use a trial wavefunction similar 
to Eq. 0, but including all four vibrational modes. The 
size of the matrix grows as 2 x C% +i m (N + 2.5) 4 /12. 
Convergence was reached for N = 10. The lowest eigen- 
values as a function of parameter a are shown on Fig. 
The coupling constants are proportional to a and at 
a = 1 they become K\ = 1.18 and k 2 = 0.8 . The fre- 
quency shift of the ring bending mode uj\ — > 0.55a->i is 
essentially the same as in the noninteracting case. The 
interaction of the bond stretching mode with the ring 
bending overtones results in the level repulsion, such that 
uj-2 — * 0.84o;2 for k 2 = 0.8. The predicted vibrational 
spectrum corresponds to a = 1 on Fig. [21 

On the anion, a D 2 h distortion splits the e 2u LUMO 
into two b\ u and a u orbitals, which can be obtained from 
Eq. (J2J by changing the sign in front of the |2 >, |4 >, 
|6 > molecular orbitals. Projecting the Hamiltonian Eq. 
Q into the {bi M ,a u } manifold leads to the same Hamil- 
tonian as J3J with a sign change of coupling constant 
gi — > The vibrational excitation spectra on the 

anion CqHq is described by Eq. with coupling con- 
stants Gi = -Tipi/6 + T 2 g 2 /2R = 0.17 eV/A and 
G' 2 = -T351/6 - T 4 g 2 /2Ro = -5.48 eV/A. Since the 
Hq part of the projected Hamiltonian is the same as 
for the cation, the values for T1..4 are unchanged. The 
JT energy lowering on the anion is the same as for the 
cation Gf/Mwl + G%/Mul = 0.234 eV, but the effective 
couplings to the ring bending and C-C bond stretching 
modes are now different; k[ — 0.16 and k' 2 = 1.06 re- 
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spectively. This means that the frequency of the ring 
bending E 2g mode is shifted by less then 2% from the 
neutral benzene value, whereas the bond-stretching E 2g 
mode becomes 59% of the neutral molecule value 1644 
cm -1 . 



VI. SUMMARY 

Using the model Hamiltonian Eq. with DFT- 

calculated parameters, we predict the JT stabilization en- 
ergy of 944 cm -1 for both the benzene cation and anion. 
The bond-stretching mode 1644 cm -1 on neutral benzene 
experiences a large shift on the anion and becomes 970 
cm -1 , which is very close in energy to the a\ g breath- 
ing mode 992cm -1 . On the other hand, the ring-bending 
mode is essentially unaffected. To the contrary, on the 
cation both E 2g JT active modes experience the sizable 
shifts by 45% and 16% from their neutral benzene values. 
Experimentally [Til |20| the ring-bending mode on the 
benzene cation was measured to reduce from 536cm -1 to 
350cm -1 (or by 35%). This measurements led Linder et 
al. Q to deduce a JT energy stabilization of 208 cm -1 . 
A more detailed analysis, including quadratic JT cou- 
pling, led to the conclusion that the acute D 2 h distorted 
geometry is the global minimum, lower by 8 cm -1 than 
the obtuse local minimum. Our Hamiltonian |0J predicts 
a larger relative shift of the ring bending mode by 10% 
with respect to the experimental shift. The amount of 
the shift is sensitive to the choice of normal modes deter- 
mined by the term Ho (Eq. 0J, a slightly oversimplified 
model. The main result of our new solution is that there 
are two contributions to the JT energy lowering. To de- 
termine the JT energy lowering on the cation two vibra- 
tional E 2g frequencies have to be measured. The model 
Hamiltonian Q contains quadratic couplings, which pre- 
dict an acute global minimum, while the accuracy of the 
true DFT answer is insufficient to make such a prediction 
□ 



The benzene anion is not stable, which makes it diffi- 
cult to measure the vibrational spectrum. Electron trans- 
mission techniques |2l| give an opportunity to measure 
the vibrational sidebands in the resonant electron cross 
section at negative electron affinity energy 1.1 eV plj . 
The spectrum shows vibrational peaks separated by 123 
meV (or 992 cm -1 ), and were attributed to a\ g vibra- 
tional quanta. The higher frequency JT-active e 2g mode 
has not been resolved. Our results, however, suggest that 
since both breathing and C-C bond stretching modes are 
coupled to the LUMO, a mixture of the two is expected 
to be present in the spectrum. Since the e 2g modes be- 
comes very close in energy to the breathing a\ g mode, 
it is difficult to distinguish experimentally between the 
two. 

In conclusion we derived a new model Hamiltonian for 
the JT active benzene cation and anion, and obtained all 
parameters from DFT calculations. This model Hamil- 
tonian predicts a JT stabilization energy of 970 cm -1 . 
This value, similar to other first-principles results, is 4.5 
times larger than the value deduced experimentally by 
Linder et al. 14]. However, our model gives the re- 
sult that two E 2g modes are JT-active, coupling to the 
electronic HOMO and LUMO states. Only one of them 
was measured and our results are reasonably con- 
sistent with that measurement. However, both modes 
are needed to obtain the JT stabilization energy from 
spectroscopic data. We conclude that there is no contra- 
diction between theory and experiment, and predict the 
relevant frequency shifts that can be used to test theory 
more completely. 
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